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ABSTRACT 

The principles of measuring the shapes of galaxies by a model-fitting approach are discussed 
in the context of shape-measurement for surveys of weak gravitational lensing. It is argued that 
such an approach should be optimal, allowing measurement with maximal signal-to-noise, 
coupled with estimation of measurement errors. The distinction between likelihood-based and 
Bayesian methods is discussed. Systematic biases in the Bayesian method may be evaluated as 
part of the fitting process, and overall such an approach should yield unbiased shear estimation 
without requiring external calibration from simulations. The principal disadvantage of model- 
fitting for large surveys is the computational time required, but here an algorithm is presented 
that enables large surveys to be analysed in feasible computation times. The method and 
algorithm is tested on simulated galaxies from the Shear TEsting Program (STEP). 

Key words: Gravitational lensing - methods: data analysis - methods: statistical - techniques: 
miscellaneous 



1 INTRODUCTION 

Measurement of the effects of weak gravitational lensing has be- 
come a key technique in the arsenal of methods used to measure 
the distribution of matter, both associated with individual objects 
such as galaxy clusters or individual galaxies, and on large-scales 
through the measurement of 'cosmic shear'. A key advantage of 
such measurement is that it directly measures the total matter distri- 
bution, generally dominated by the dark matter component, which 
may then be related directly to theory without needing to under- 
stand the uncertain effects of the physics of baryons in galaxies , 
provided one avoids the highly nonlinear regime dWhitel 12004 . 
IZhan&Knoxll2004 Jing 2006). Thro ugh the use of photometric 
redsh i fts, three-dimen sional analyses jHul 1 19991 : iBacon & TavloJ 
120031 : lHeavensll2003h can be used to further measure both the 
cosmologic al growth of structure and the values of cosmological 
parameters (Massey et al. 2007a; Heavens et al. 2006; Taylor et al.l 
l2007l : lKitching et alj|2007h . Until recently such surveys have been 
of limited size, but even so the results obtained provided useful 
constraints on cosmological parameters and an important test of 
the values deduced from other methods. One long-standing puzzle 
has been that the range of values for the power-spectrum normal- 
isation parameter as found by weak lensing analyses has tended 
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to be higher than found b y some other methods (see the discus- 
sion in lSpergel et al.ll2007h . an effect that persists at some level in 
the latest studies. For the best fit 3-year WMAP value of the mat- 
ter dens ity parameter Qo = 0.24, the 3D analysis of Mass ey et al.l 
(2007b |) finds the valu e a s = 0.96±;°7, and the 2D analysis of 
Benjamin et al. (2007) finds as = 0.84 ± 0.07. These results can 
be compared with the 3-year WMAP value as = 0.76 ± 0.05 
dSpergel et alj2007h . 

Measurement of the effect of weak gravitational lensing re- 
quires the statistical analysis of large samples and is sensitive 
to any systematic errors in measured quantities. Possible system- 
atic errors in lensing signals introduc ed by uncertainty in photo- 
metric redshifts has been discussed by Edmonds on. Miller & Wolj 
(2006). Another fundamental concern with the method is whether 
the shapes of galaxies, that are used to deduce the signal, may be 
measured in an unbiased manner. The problem of shape measure- 
ment in optical imaging data is that galaxy images are convolved 
with a possibly-varying point-spread function (PSF) which must 
be accurately corrected for when deducing galaxy shape. Convolu- 
tion with the PSF tends to make galaxy images appear rounder (for 
reasonably circularly symmetric PSFs) whereas addition of pho- 
ton shot noise has the systematic effect of tending to make round 
galaxies appear less round. These two observational effects thus 
tend to work in opposite senses, and are independent of each other, 
so that both accurate PSF correction and calibration to remove the 
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effects of noise on shape are required. Fol lowing the seminal pa- 
per bv lKaiser. Squires & Broadhurstl d 19951) there have been many 
su ggestions for po s sible m eas urement proce s ses, tha t are discussed 
bv lHevmans et al.l d2006l) and lMassev et alj d2007bl) as part of the 
'Shear TEsting Program' (STEP). Those papers discuss 18 pub- 
lished methods for shear measurement. The existence of so many 
suggested methods implies that no consensus has yet emerged on 
the best way to measure weak lensing signals, and therefore nat- 
urally leads us to ask whether there might in fact be one method 
that may be regarded as being optimal. In this paper we investigate 
whether a model-fitting approach to galaxy shape measurement can 
both achieve this aim of optimal measurement and also be con- 
structed such that it is computationally feasible for large surveys. 

In the following we shall suppose that galaxies may be char- 
acterised by a measurement of their ellipticity e and that a weak 
lensing signal, such as cosmic shear, may be inferred either from 
the mean ellipticity or from some form of cross-correlation of the 
ellipticities of different galaxies. 

We can stipulate a number of requirements that a weak lensing 
measurement technique should satisfy. 

(i) Optimal measurement of lensing signal, in the sense of max- 
imum signal-to-noise. 

(ii) Unbiased measurement of lensing signal. 

(iii) Ability to calculate the statistical uncertainties of the mea- 
surement. 

A standard approach that in principle allows us to meet these crite- 
ria is that of model-fitting, which is the method discussed in this pa- 
per. We first discuss some general principles, including whether we 
should use a frequentist or Bayesian approach and how shear may 
be measured in an unbiased way from a Bayesian posterior proba- 
bility distribution. However, the principal disadvantage of a model- 
fitting approach is that it might be computationally prohibitive for 
very large surveys. In section[3]we discuss a novel galaxy shape 
model fitting algorithm that allows good estimation of the likeli- 
hood surface in a usefully short computational time. We also dis- 
cuss the evaluation of shear sensitivity within the Bayesian frame- 
work, that allows individual galaxy contributions to be assessed and 
unbiased estimation of shear to be made, fulfilling the second crite- 
rion above. Some initial results and further considerations are then 
discussed. More detailed results from applying the algorithm to the 
STEP simulations are given in a companion paper (Kitching et al., 
in preparation). 



2 A MODEL-FITTING APPROACH TO SHAPE 
MEASUREMENT 

2.1 General considerations 

The basic rationale for fitting a model of a galaxy's surface bright- 
ness distribution is that, if the family of models is a good repre- 
sentation of the true surface brightness profile, the highest possi- 
ble signal-to-noise of the resulting parameters should be obtained. 
When model and data agree the model encapsulates the full in- 
formation content of the data. Although this h as been recognised 
previ ously in weak lensing shape measurement 1 Berns tein & Jarvisl 
2002), no implementation of weak lensing shape measurement 
methods published to date has this property, because the meth- 
ods usually adopt some simplification of the surface brightness 
profile, such as assumi ng that second momen t s entirely charac- 
terise the profile (e.g. iTvson. Wenk & Val des 1990 Kais eretai] 



1 19951) or equiva l ently assuming Gau s sian profiles or w eights 
(e.g. iBridle et al] 120021 : iKuiikenl Il999l : iBernstein & Jarvisl |2002|; 
Bard eauetal .1120051) . Model-fitting has been used for some time 
for detailed determination of g a laxy surfa ce brightness profiles and 
shapes (e.g. IPeng et al.l 1200 3). Kuijkerj d 1999b proposed model- 
fitting to avera ged galaxy images specifically for weak lensing mea- 
surement, and lBridle et al.l 112002) proposed a method of measuring 
shear by fitting galaxies and PSFs with multiple Gaussian compo- 
nents. The latter method has been applied to survey s of weak lens- 
ing aro und galaxy cluste r s by Bardeau et al.l d2005l) . iBardeau et al.l 
(2007) and iKneib et al] d2003l) , among others. A Monte-Carlo 
method is used to find best-fitting galaxy model parameters for each 
individual galaxy, where Gaussian surface brightness profiles, or 
combinations of two Gaussian profiles, are assumed for both galaxy 
and PSF. Shear measurement and the computation al time required 
for tha t model-fitting method has been evaluated bv lHevmans et al] 
(2006). Recently, sets of basis functions known as 'shapelets' have 
been used to describe su rface brightness profil es |Refregier]|2003l: 
iRefregier & Bacorj2003l and the related work of Bernstein & Jarvisl 
2002) but there is no requirement for the individual shapelet func- 
tions to match real galaxy profiles. Moving to a pure model-fitting 
approach allows us to choose whichever brightness profiles we like, 
and for galaxies it clearly makes most sense to choose either ex- 
ponential or de Vaucouleurs surface brightness profiles. Naturally, 
the above statements are qualitative, we don't know how much im- 
provement one obtains by fitting a profile that is closer to the actual 
profile, but the principle at least is a sound one, that we expect to 
satisfy the first of our criteria from the Introduction. 

Either frequentist model-fitting, based on determining the 
likelihood function, or Bayesian model-fitting that determines the 
posterior probability distribution of model parameters, allow error 
estimates to be made, satisfying the third of our criteria. This is 
not the case for early versions of weak lensing shear estimators, 
althou gh error estimate s have been made in some recent meth- 
ods dBernstein & Jarvisl [ 20021 : IKuiikenl I2OO6I IBridle et al]|2002l: 
IBardeau et al]|2005l) . 

Finally, we should address the question of whether a method 
can be determined to be unbiased. This is a serious issue for weak 
lensing studies, where the signal is so small that even a small sys- 
tematic bias can have a devastating effect. Evaluation of existing 
methods by STEP demonstrate that they are indeed biased, with 
significant magnitude-dependent biases that need to be corrected 
empirically from compar ison with simulations dHevmans et al] 
2006 ; iMassev et alj |2007b ) . We discuss in the next section why in 
principle a Bayesian method should be unbiased provided a cor- 
rect choice of prior is made, but note that realistic implementations 
result in a quantifiable bias that may be corrected for. 



2.2 Bayesian estimation of the sample ellipticity distribution 

We have previously mentioned the problem of shape measurement, 
that not only is the shape changed by convolution with the PSF, but 
also noise biases the measured shape, and in general tends to make 
nearly-circular objects systematically appear more elliptical. We 
discuss in this section how a Bayesian method may be formulated 
that precisely corrects for this phenomenon, provided we make a 
correct choice of prior. 

Consider a set of observations of N galaxies that yields the 
surface brightness distribution for each galaxy denoted by a vector 
of pixel values y . The shape of each galaxy may be characterised by 
its two-component ellipticity e: the particular definition we choose 
for e in this paper is given in section l2~4l but what follows below 
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applies to any shape estimator that we may choose. If the sample of 
galaxies has a probability distribution of intrinsic ellipticities (i.e. 
the value of ellipticity that would be measured by the observer in 
the absence of degradation by the PSF or by noise) /(e), then the 
probability distribution of y is 



n(y) = J f{e)e{y\e)de, 



where e(y e) is the probability distribution for y given e. 

For each of these galaxies we can generate a Bayesian poste- 
rior probability distribution for its ellipticity 



Pi(e\Vi) 



P(e)C(y 1 \e) 
JP(e')C( yi \e')de' 



where V (e) is the ellipticity prior probability distribution and 
£ (yje) is the likelihood of obtaining the i th set of data values 
y i given an intrinsic ellipticity e. 

We would hope that the true distribution of intrinsic elliptic- 
ities can be obtained from the data by considering the summation 
over the data: 



i 

V(e)£(y\e) 



I 



Ay 



JV(e')L(y\e')de' 



f(e")e(y\e")de' 



where on the RHS we are integrating over the probability distri- 
butions to obtain the expectation value of the summed posterior 
probability distribution for the sample. We can see that this will be 
achieved if both e(i/|e) = C (y\e) and V (e) = / (e), assuming 
the likelihood is normalised, J C (y\e) dy = 1, from which we 
obtain 

i 

The strength of this result is that we can in principle recover 
statistically knowledge of the intrinsic distribution of shapes inde- 
pendently of assumptions about the shapes of the likelihood sur- 
faces: in particular the likelihood surfaces for ellipticity measure- 
ment must be non-Gaussian, being bounded at |e| < 1, but this 
has no effect on the results we expect. This result parallels the 
analogous result discussed by Edmondson et al. (2006) for the case 
of Bayesian photometric redshift estimation. It says that we must 
know the mechanism by which data values are generated in or- 
der to construct the likelihood function, and that we must know 
the expected distribution of intrinsic ellipticities, in which case the 
summed posterior probability distribution will recover that intrinsic 
distribution. It might be thought that a Bayesian approach then has 
a non-useful requirement, that we need to know the answer before 
we start, but the point of course is that with the correct choice of 
prior we then expect the posterior probability distribution for each 
individual galaxy to yield an unbiased estimate of ellipticity, and 
those sets of individual posterior probability distributions may then 
be used to infer the spatially varying shear arising from gravita- 
tional lensing. We discuss in section |5T| one possible method for 
creating the correct prior. 



2.3 Frequentist or Bayesian measurement? 

So far the framework has been described in a purely Bayesian con- 
text, but we can also ask whether there is a frequentist equivalent 
of the above formalism: can weak lensing shear be measured using 
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Figure 1. Illustration of the properties of an ideal likelihood estimator xc 
(left) and ideal Bayesian estimator xg (right) for the Gaussian example de- 
scribed in the text. The top pair of graphs show the correlation between 
the input and deduced values. Two regression lines are shown on each, 
one being the regression of input on estimated value, the other being the 
regression of estimated on input value. The next two pairs show the dis- 
tribution of the diff erence between input and estimated values compared 
with either the estimated values (centre) or the input values (bottom). Note 
the graph x-axes differ between the centre and bottom panels. For a given 
input value, the likelihood estimator yields an unbiased estimate (regres- 
sion slope unity) whereas the Bayesian estimator appears biased (regression 
slope 2.75). However, for a given estimated value, the likelihood estimator 
is biased (regression slope 0.36) and the Bayesian estimator is unbiased. 
The Bayesian estimator returns the best estimate of the input value for a 
given measurement. 



likelihood functions alone? Conversely, are there any disadvantages 
to using a Bayesian method? 

It is important to recognise that likelihood and Bayesian esti- 
mators measure different things. We can illustrate this by consider- 
ing a sample of galaxies (say) with some intrinsic property x that 
we wish to determine from fitting to some measurements y. We 
shall look at the results obtained with either a Bayesian estimator, 
xb = / xp(x | y)dx, or a likelihood estimator, xc = J xC(y\x)dx 
(the general considerations discussed here apply also to maximum 
likelihood estimators). Suppose the intrinsic distribution of x has 
a normal distribution of variance a 2 , and that for each x drawn 
from this distribution the measurement process causes a normally 
distributed uncertainty of variance b 2 . Fig.Q] shows the results ob- 
tained in a Monte-Carlo realisation for the illustrative case a — 0.3, 
b = 0.4. 



© 2007 RAS, MNRAS 000. [71071 



4 L. Miller et al. 



The likelihood estimator is based on the function £(y\x) and 
in Fig.Q]is unbiased in the regression of input on estimated val- 
ues. The Bayesian estimator is based on the function p(x\y) and 
is unbiased in the regression of estimated on input values: that is, 
for a given set of measurements the Bayesian estimator yields the 
best estimate of the input values. The likelihood estimator yields 
a distribution of measured values that is broader than the intrinsic 
distribution. Notice however that the Bayesian estimator yields a 
distribution that is narrower, despite the result of section [2~2l that 
the summed posterior probability distribution yields the intrinsic 
distribution if the correct prior is chosen. This apparent paradox is 
resolved by realising that each of the individual estimated values is 
associated with its own posterior probability distribution, so that the 
sum of the distributions is broader than the distribution of expecta- 
tion values. While this might seem undesirable, this is inevitable in 
any noisy measurement process. 

In the case of ellipticity shape measurement, we expect there 
to be other reasons why a likelihood estimator might be biased, in 
particular at large e values or for low signal-to-noise, where the 
boundary |e| < 1 renders the likelihood function asymmetric and 
highly non-Gaussian. With no prior, the strong degeneracy between 
size and ellipticity in the likelihood fitting can create regions of 
high likelihood at extreme values of ellipticity, and given the hard 
bound |e| < 1 such an estimator cannot be unbiased. 

In the frequentist approach, it is also possible to estimate er- 
rors for the individual galaxies, and this is important to establish 
the contribution to the signal from each galaxy. As in the Bayesian 
approach, we expect the information on shear to decrease as the 
signal-to-noise (S/N) decreases: even with an unbiased estimator e 
it would be important to quantify this effect and allow for it in the 
shear estimation. 

A final consideration is that in weak lensing surveys we are 
not simply interested in measuring the shapes of individual galax- 
ies, but rather in measuring the systematic lensing shear in a sam- 
ple. Usually this is measured from the mean ellipticity: thus it may 
be possible to have an ellipticity estimator that is biased but where 
the shear estimate from a sample is unbiased, or vice versa. If any 
bias were isotropic, corresponding to a bias in the value of e but 
not in orientation, then we might hope that the bias would aver- 
age out. However, even in this case we should not assume that the 
shear estimator (e) is unbiased, since the likelihood functions for e 
measurement must be non-Gaussian and e-dependent, any shift g 
in the distribution of e caused by lensing would lead to bias in the 
estimated shear. 

However, even the Bayesian method is not immune to the 
problem of bias, particularly in a realistic implementation of the 
Bayesian method where we are forced to assume a zero-shear prior, 
as discussed below. But the bias can be quantified and the method 
provides a self-contained framework within which we can work out 
all the required quantities. This is the framework that we return to 
in the remainder of this paper. 



2.4 Bayesian shear estimation and the shear sensitivity 

Following Hevmans et al. ( 2006) we assume observed galaxy ellip- 
ticity e is related to the intrinsic galaxy ellipticity e 3 in the weak 
lensing regime via: 

e 3 +g 
e i+g*e 3 

from lSchramm & Kavsej dl995t) : ISeitz & Schneidel ( Il997» . where 
e is represented as a complex variable and g, g* are the reduced 



shear and its complex conjugate respectively, e is defined in terms 
of the major and minor axes and orientation a, b, 6 respectively, as 
e = (a — b)/(a + b) exp(2i#). In this formalism, we expect 



(1) 



for an unbiased sample where (e s ) = 0, and so (e) for a sample 
of galaxies is adopted as our estimator of g. Note that this result 
differs from the other commonly used formalism where ellipticity 
is instead denned as e = (a 2 — b 2 )/(a 2 + b 2 ) exp(2i#). 

For a population of galaxies, (e) = J e/(e)de where /(e) 
is the ellipticity probability distribution for the sample. But in the 
Bayesian formalism we can write a similar expression for an indi- 
vidual galaxy if we know its Bayesian posterior probability distri- 
bution, {e)i — J ep(e|y i )de. Hence for a sample of TV galaxies 
we can evaluate the sample mean as 



In practice we shall use the first of these expressions, as estima- 
tion of ellipticities for individual galaxies allows error estimates to 
be made for each galaxy, and its contribution to the signal to be 
evaluated. 

However, in measuring shear we cannot know in advance the 
correct prior to apply, even if we know the intrinsic unsheared el- 
lipticity prior distribution, because the amount of shear varies over 
the sky in a way that we are attempting to measure. We must there- 
fore use a prior that contains zero shear. The effect of this is that 
as signal-to-noise decreases, the measured ellipticity distribution 
tends to the prior, and in the limit of zero signal-to-noise no shear 
signal is recoverable. This is precisely what we should expect of 
course: no measurement method can extract a measured shear value 
from data with zero signal-to-noise, and a Bayesian method is no 
different in that respect. A Bayesian method does however allow us 
to estimate the magnitude of this effect for each individual galaxy. 
Consider the Bayesian estimate of ellipticity (e)i, denned above, 
that is measured for the i th galaxy, and express its dependence on 
each component of shear g as a Taylor series. For component e\, 



(ei)i ~ e\i + gid{ei)i/dgi + g 2 d{e 1 ) l /dg2 + ■ 



(2) 



and similarly for component e 2 , where numeric subscripts indicate 
the components of e and g. In the weak lensing limit the cross- 
terms vanish. If we sum over TV galaxies in an unbiased sample we 
find 



We may optionally multiply both sides in equationf2]by a statistical 
weight for each galaxy, Wi . Provided Wi is uncorrelated with ef we 
may then define a weighted estimate of shear for the sample: 



(3) 



E, ^d{e M )»/<9g M 

for fi — 1,2. We shall call d(e l± )i/dg ti the shear sensitivity. It 
is a measure of how much each Bayesian estimate is biased by 
the use of the zero-shear prior, and it takes values in the range 
< d{e )J )i/dg li ^ 1, where the lower bound is expected in the 
limit of zero signal-to-noise. The upper bound would be attained 
in the case of ideal measurement at high signal-to-noise, where we 
expect no bias: in this case the Bayesian measure is a good measure 
of the true ellipticity, regardless of which prior is assumed, and dif- 
ferentiating equationQ] yields unity for the shear sensitivity. The 
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weights Wi could in principle be tuned for optimal signal-to-noise 
in the measurement(s) being made, such as the values of cosmo- 
logical parameters. Care should be taken that any weights that are 
a function of ellipticity do not introduce bias into shear measure- 
ment. Since the shear sensitivity is effectively a measure of how 
much information about the effect of lensing is carried by each 
galaxy, the weights should themselves also include a dependence 
on d{e l _ l )i/dg fl , as well as on the measurement error and on the 
redshift-dependent cosmological effect of lensing on each galaxy. 

The shear sensitivity may be estimated for each galaxy and 
for the survey as a whole without recourse to external calibra - 
tion from simulat i ons, a s describ e d bel ow. iKaiser et alj dl995l): 
iLuppino & Kaiserl dl997h : I Raised feOOd) and iBernstein & Jarvij 
l 2002h have emphasised the utility of knowing the shear 'polar- 
izability' or 'responsivity' for individual galaxies, as this not only 
allows accurate optimised shear measurement but also allows fu- 
ture surveys to be planned and optimised. 

The estimator g^ is appropriate for a survey where the shear 
is uniform over some region (and this is assumed in the STEP sim- 
ulations discussed below), but in the more general case we instead 
infer the shear correlation function or some related quantity such as 
shear variance from a measurement such as (ejej). In this case we 
can compute the analogous estimator 



3^) = 



We now discuss possible approaches to calculating the shear 
sensitivity, first for normal prior and likelihood distributions, then 
for the more general case where the shear sensitivity may be evalu- 
ated numerically from the measured likelihood surfaces of individ- 
ual galaxies. 

2.5 Calculation of shear sensitivity 

As an illustration of the calculation of shear sensitivity, suppose 
the prior is described by a normal distribution V(e) of variance 
a 2 and (e) = 0, and that the likelihood £(e) for a particular 
galaxy also has a normal distribution of variance b 2 centred on 
some value eo. It is straightforward then to show that the Bayesian 
posterior probability p(e\y) also has a normal distribution of vari- 
ance a 2 b 2 /(a 2 + b 2 ) and expectation value (e) = eoa 2 /(a 2 + b 2 ). 
For perfect measurement of ellipticity (b 2 <C a 2 ) we expect equa- 
tion[T]to hold, so for this galaxy d(e li )i/dg fl = deo^/dg^ = 1. 
For more noisy measurement, we expect 



<9{e p 



<9e . 



%m « 2 + b 2 dg^ a 2 +b 2 ' 
The shear sensitivity decreases as the measurement error increases. 
The value of the shear sensitivity is also given by the inverse of 
the slope of the regression of the intrinsic ellipticity on estimated 
ellipticity illustrated in Fig.[T] 

The above illustration indicates that it is straightforward to 
calculate the shear sensitivity, however in general it would not be 
safe to assume normal distributions: not least because e is defined 
such that |e| < 1, so when the measurement error becomes large 
£(e) cannot be normally distributed. We discuss here one method 
of calculating the shear sensitivity numerically. We should empha- 
sise that this can be done entirely internally to the fitting process, 
with no need to calibrate shear sensitivity externally from simula- 
tions. 

Consider first the response of the posterior probability distri- 
bution to a small amount of shear. The prior probability does not 



depend on the shear in our implementation. Let us assume that ap- 
plying a weak lensing shear shifts the likelihood function by some 
small amount, £(e — e 3 ) — > Lie. — e" — g) and expand as a Taylor 



£(e - e" - g^ 
Then, substituting into 



9n 



dC 

de u 



JeV{e)£{e)de 



fV(e)£(e)de 
and differentiating with respect to g we find, 

<9(e M 



f((e)-e)V(e)§±de 



JV(e)C(e)de 



as an estimate of weak lensing shear sensitivity. This expression is 
cast in terms of the derivatives of the likelihood surface multiplied 
by the prior: it may also be expressed in terms of derivatives of the 
prior multiplied by the likelihood: 



<9{e p 



dg^ 



1 



J((e)-e)C(e)^de 
JV(e)C(e)de 



This may be evaluated numerically from the posterior probability 
surface for each galaxy, and is preferred over the preceding expres- 
sion in the case where the derivative of the prior is known analyt- 
ically. For the case of normal distributions of V(e) and £(e) the 
expression yields the analytic result above. 



3 FAST SHAPE MEASUREMENT 
3.1 The algorithm 

The technique we adopt for measuring (e) and its uncertainty is to 
fit model galaxy surface brightness profiles to the data for individ- 
ual galaxy images. The simplest galaxy model has six free param- 
eters if the form of the surface brightness profile is fixed: central 
surface brightness, size, ellipticity and celestial position. The prob- 
lem of fitting six parameters to large samples of galaxies is that this 
could be a time-consuming task, probably prohibitively so. How- 
ever, we can greatly speed up the process if we can marginalise 
over any parameters that are not of interest to the weak lensing mea- 
surement. It turns out that for isolated galaxies it is straightforward 
to marginalise over three of the parameters, central surface bright- 
ness and position, if the model fitting is treated in Fourier space, 
as described below. And because there exist fast Fourier transform 
algorithms this approach can be done in a short amount of compu- 
tational time. 

We can start by writing the statistic 



E 

i 

£ 



A(C - B) 2 - AB 2 , 



where yi is the data value in pixel i, Oi is the statistical uncertainty 
of that data value, y™ is a model value for that pixel, C is the model 
amplitude and where 

Vi \ B _^ViVi 

Oi 



A 



£ 



E^h/E 



We assume the pixel noise is stationary and uncorrelated, which is 
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appropriate for shot noise in CCD detectors in the sky-noise limit. 
Bright galaxies may make a significant contribution to photon shot 
noise, but in this case the non-stationarity of the noise makes it not 
possible to work in Fourier space. Hence this algorithm is appropri- 
ate for model fitting to faint galaxies in the sky-limited regime. The 
method can be generalised to the case where the noise is stationary 
but correlated between pixels, and the example of radio interferom- 
eter observations is discussed qualitatively in section [5l4l 

Then if we adopt a prior V(C) for the model amplitude we 
can marginalise the likelihood L = e~ x I 2 over C: 



C = e 



-T,vt/ 2 ^ e AB / 2 



-A(C-Bf/2 nc)da 



We shall adopt a uniform prior for V(C) in the range C m m ^ C ^ 
Cmax- We expect C > 0, but if the galaxy is significantly detected 
the Gaussian form of the likelihood causes the value of the prior to 
become unimportant at both large and small C, and we can simplify 
the calculation by allowing C m i n — > — oo and C" m ax —> oo, so that 



However, although we have eliminated the amplitude C, C still de- 
pends on the model second moment, A — E (y™/<7i) 2 . Thus we 
need to introduce the model constraint A ^constant, achieved by 
renormalising each model appropriately. Since for a given dataset 
^2yf /of is also fixed, we can write C oc e AB I 2 when maximis- 
ing. 

We can also rapidly calculate the marginalisation over galaxy 
celestial position if we work in Fourier space, writing 



E 



J/fce 



y, 



E 



Vk e 



We can simplify the various summations by assuming that we are 
dealing with faint galaxies in weak lensing measurement, such that 
<7i is dominated by the background photon shot noise and is con- 
stant for all pixels. And since the model y™ is real, = y™* and 
2~2i ViVT — 2~2k VkVk 1 * ■ If we introduce a shift X into the model 
position, the new model becomes 

m/ \ m —ik.x; —ik.X 

Vi = z^Vk e e 



and 



ViVi =2^ykVk e = h(X) 



where h(X) is the cross-correlation of the data yi with the model 
y™- So the likelihood becomes 



C oc exp 



MX) 



2cr 2 E vT 



To marginalise over X we need to adopt a prior V(X), but in this 
case it cannot be uniform as C — >constant as \X\ — > oo and the 
marginalised likelihood would not be finite. This problem arises 
because, no matter how large a pixel value, it always has a finite 
chance of being due to random noise, with the true galaxy being 
positioned elsewhere. We shall adopt a prior which is centred on 
some assumed galaxy position that has been previously estimated 
and which falls off to zero at large distances: this is equivalent to 
assuming that a galaxy does indeed exist somewhere near the loca- 
tion we have chosen. We shall assume a prior which is symmetric 
and centred on the nominal galaxy position, which for convenience 



is at the coordinate origin, such as: 

V{X)d 2 X = -^—e-^' 2b2 d 2 X. 

The process of model fitting is seen from the above to be one 
of cross-correlating the data with a model. Galaxies generally 
have smooth centrally-concentrated surface brightness distributions 
which are convolved with near-Gaussian PSFs in an observed im- 
age. The model is also smooth, centrally concentrated and con- 
volved with the same PSF. From the central limit theorem such a 
cross-correlation should be well represented by a two-dimensional 
Gaussian distribution, 

h(X) = h exp [-(X - X )C'\X - X ) T ] 

where C' 1 is the inverse covariance matrix and some shift Xa 
of the maximum from the origin is allowed. In what follows we 
assume circular symmetry for simplicity, although this assump- 
tion may be removed without affecting the final result (a two- 
dimensional Gaussian distribution may always be transformed to a 
circularly-symmetric distribution by a coordinate transformation). 
If the cross-correlation function has the form h — ho exp[— \X — 
Xo\ 2 /s 2 ] then 



C oc 



1 



2TTb 2 



exp 



where 



ix-x r/s 2 



hi 



-\X\ 2 /2b 2 



^ 2 T.y. 



We could evaluate this by, for example, expanding the first expo- 
nential as a Taylor series and hence obtaining a series solution for 
the marginalised likelihood. We could also evaluate it purely nu- 
merically, but this would require evaluation of the cross-correlation 
function on an extremely fine grid in order to achieve adequate ac- 
curacy. Either of these approaches would be computationally ex- 
pensive, and an alternative is to find an approximate value of the 
integral by writing 



C oc 



1 



2-wb 2 
If b > s, 

£ oc 1+ 



2-wb 2 



I exp ^(3e 



-\x r/2b- 



-\X-X \ 2 /s 



]-' + '}< 



-\X\ 2 /2b' 



d z X. 



exp \(3e 



\X-X \ 2 /s' 



L|d 2 X 



and changing variables to a polar system centred on X , 

{exp [f3e- r2 ] - l} 

(3 > 1. 



1 + 



\X \ 2 /2b 2 



'/2b 2 



rdr 



213b 2 



(4) 



approximately, where the constant of proportionality has no model 
dependency provided A is held invariant (we could obtain a simi- 
lar result more exactly if we were to adopt a top-hat prior for the 
galaxy position). If the width s, amplitude ho and centroid Xo of 
the cross-correlation function can be measured, the marginalised 
likelihood may be estimated from equation|4] In the more gen- 
eral case, where the cross-correlation function is approximated by 



It may seem that the requirement for a prior on position may be removed 
by allowing b — » oo. However, this is an artefact of the approximation. 
There is no clear way of identifying a value for b, but it should be set suffi- 
ciently small that confusion from other nearby galaxies is eliminated. 
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a bivariate Gaussian with widths si and S2 in the two principal di- 
rections, then equation[4]is modified by s 2 —* S1S2. 

3.2 Implementation 

Using the above, the likelihood may be estimated for a given set 
of model parameters, marginalised over the 'uninteresting' position 
and brightness of the galaxy. If we choose either an exponential disc 
or a de Vaucouleurs model for the surface brightness, the free pa- 
rameters are the scale-length and two ellipticity values e, or equiv- 
alently the scale-length, the axial ratio and the orientation. The val- 
ues of e are restricted to lie in the range |e| < 1, and for faint 
galaxies the probability p(e) is broad, making a grid search in e an 
easy and not-too-expensive approach. The resulting likelihood may 
be numerically marginalised over the galaxy scale-length, which 
is also 'uninteresting' for weak lensing measurement, to obtain a 
likelihood surface that is a function of ellipticity alone. 

To evaluate the cross-correlation function h(X) (returning 
now to the general case where we do not assume circular symme- 
try) we can use the fast Fourier transform method, which proceeds 
as follows: 

(i) Generate a series of 2D galaxy surface brightness models on 
a three-dimensional grid in parameter space of scale-length and 
ellipticity. These models can be discrete Fourier transformed and 
those transforms stored for use with all the galaxies. The choice 
of a grid of models allows a considerable multiplex gain to be re- 
alised: the models can be pre-generated on that grid and the same 
set used for fitting to every galaxy. 

(ii) Estimate the surface brightness profile of the PSF on the 
same pixel scale as the models. Usually this would be done by 
stacking images of stars from the region of an image as the galax- 
ies being measured. The PSF can also be Fourier transformed and 
stored. If the PSF varies over an image or between images, the im- 
age may be divided into zones over which the PSF is approximately 
invariant, and the Fourier transform of the PSF for each zone stored 
separately. If a mathematical model for the varying PSF is known 
this may also be u sed to generate a smoothly varying PSF (e.g. 
iRhodes et al]|2007t) . 

(iii) Estimate the rms noise in each pixel from the entire image. 

(iv) Identify a set of nominal galaxy positions to be measured, 
most likely from a separ ate image analysis tool such as SExtractor 
l lBertin & Arnoutsll 19961) . 

(v) In turn for each galaxy, extract a sub-image centred on that 
galaxy, Fourier transform it, and temporarily store the result. 

(vi) For this galaxy, take each possible model in turn, multiply 
by the transposed PSF and model transforms to carry out the cross- 
correlation, measure the amplitude, width and position of the maxi- 
mum of the resulting cross-correlation, and hence evaluate the like- 
lihood for this model and galaxy. Repeat for all models on the grid 
(or for a subset of models if a more intelligent maximum-likelihood 
or MCMC search algorithm is being employed). 

(vii) Numerically marginalise over the scale-length parameter. 
In the implementation described here we adopt a uniform prior for 
the distribution of galaxy scale-length. This could be replaced by a 
prior close to the actual distribution of galaxy sizes, although such 
a prior would need to be magnitude-dependent. 

(viii) Discard the extracted data when all models have been ex- 
plored, and repeat for the next galaxy. 

The result is a grid of likelihood values in ellipticity parameter 
space which thus defines the probability surface p(e). The reduced 
shear may then be directly estimated from (e), and the uncertainty 



in individual e values may be estimated from the width of the like- 
lihood surface. 

There is a significant multiplex gain obtained by Fourier- 
transforming the models, the PSFs and the data and storing the 
results. The time-consuming step then is the cross-correlation, 
which comprises some multiplications and a single inverse Fourier 
transform to obtain the cross-correlation function. It is this multi- 
plex gain, combined with the elimination of three parameters by 
marginalisation, that yields a fast fitting algorithm. 

The algorithm is approximate, in the sense that we require the 
cross-correlation amplitude to be high enough that (3 ^> s 2 /2b 2 , 
and also in that we assume the core of the cross-correlation function 
can be adequately modelled as a Gaussian, and we have assumed 
that the pixel noise is invariant. This latter constraint may impose 
a maximum brightness limit on galaxies that may be fitted, as the 
pixel noise is not invariant in the case where the galaxy itself makes 
a significant contribution to the noise. For a fixed size of extracted 
region around each galaxy, there is also a maximum galaxy size 
that can be adequately measured. Larger sizes are possible at the 
expense of greater computation time. 

We note that, in this method, the final PSF that is used is it- 
self a convolution of PSF components arising from the atmosphere, 
telescope and the pixel response of the detector. We do not need to 
distinguish the origin of the final PSF that is used, the method takes 
a galaxy model and convolves that with an estimate of the final 
PSF in order to cross-correlate with the data. Ultimately however 
this, and all shape measurement methods, are limited by the extent 
to which the sampled data fully encapsulate the information in the 
sky: the effect of sampling is to alias spatial frequencies higher than 
the Nyquist sampling frequency. This affects both the creation of 
the stacked PSF and the model-fitting itself. If astronomical obser- 
vations were band-limited this would not be a problem, but in real- 
ity some aliasing is inevitable. Poorly sampled observations should 
ideally be "dithered" in order to reduce such aliasing effects. 



4 RESULTS 

4.1 Tests on simulated galaxy images 

The algorithm has been implemented and tested on simulations 
provided for the 'Shear TE sting Program', STEP dHevmans et al.l 
l2006l : lMassev et al]|2007bl) . Images of galaxies were simulated for 
the Canada-France-Hawaii telescope with pixel scale 0.206". The 
simulations used here to demonstrate basic shape measurement are 
those with an isotropic PSF of FWHM 0.9" and zero lensing shear 
(tests of shear measurement in the companion paper will cover all 
the simulated PSF shapes and shear values). Simulated galaxies 
with a mixture of bulge/disc components were used but all were 
fitted with a single exponential surface brightness profile. 

As here we are testing the Bayesian method, and not our abil- 
ity to locate galaxies, we use as input galaxy positions those that 
were used when making the simulations. We also adopt as the prior 
V(e) the input ellipticity distribution used in the simulations. 

For these tests the size of each subimage was 32 pixels square. 
The choice of subimage size is a compromise between (i) having 
the subimage large enough that the galaxy surface brightness distri- 
bution is not unduly truncated and (ii) not allowing the computation 
time to become excessively long. In our initial implementation we 
have also required that only a single galaxy should occupy each 
subimage, thereby eliminating close pairs. The choice of 32 pixels 
for the STEP galaxies ensured that the subimage was larger than 
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Figure 2. Comparison of Bayesian posterior probability p(e) (left) and 
likelihood C(e) (right) surfaces for two individual galaxies. The grey-scale 
is logarithmic showing a range of 5 in A log C below the maximum value 
(shown as white) in each case. The upper panel shows results from fitting 
to a magnitude 24.17 simulated STEP galaxy, the lower panel a magnitude 
23.15 galaxy. Solid lines show the two parameter 1-cr and 2-<r contours. 
The cross shows the input ellipticity value. 



the half-light diameter in every case. In principle, the subimage size 
could be a function of galaxy size or brightness, but this sophisti- 
cation would introduce some complexity into the code and has not 
been tested here. 

The likelihood was evaluated on a cartesian grid in e, sam- 
pling at intervals of 0.1 out to a maximum axial ratio a/b = 10. 
The resulting galaxy shapes were found to be consistent for elliptic- 
ity grid intervals less than 0. 1 for these STEP galaxies. The choice 
of grid interval may need to be adjusted for different surveys. 

The PSF was created by stacking stars from the simulation, al- 
lowing sub-pixel registration using sinc-function interpolation. Ul- 
timately any shear measurement survey will be limited by the ac- 
curacy to which the PSF is known. Systematic PSF errors will of 
course cause a systematic error in estimated shear, and if the PSF 
varies on some angular scale within a survey this will imprint a 
signal on that scale on the shear power spectrum. This concern is 
common to all methods of shape and shear measurement, and we 
do not specifically address this problem here. 

An assumption of the fast fitting algorithm is that we are fitting 
to individual galaxies, and hence close pairs of galaxies cannot be 
fitted with this algorithm. In practice one could identify such close 
pairs in the data at the galaxy-detection stage, and on those galaxies 
we could use a fitting algorithm that fits multiple components. In 
this case the full six parameters per galaxy would need to be fitted, 
with marginalisation over uninteresting parameters being carried 
out post-fitting. There would still be a significant speed advantage 
to be gained by using the fast fitting algorithm on the more isolated 
galaxies however. In this paper we focus on testing the Bayesian 
method and the fast fitting algorithm, and hence in the results pre- 
sented here we exclude cases where multiple objects are identi- 
fied within a single galaxy sub-image. This procedure excludes 



Figure 3. Tests on the STEP 1 simulated galaxy sample, as a function of 
galaxy apparent magnitude. Each graph shows the expectation value of the 
Bayesian estimate of component e\ (x-axis) plotted against the input value 
(y-axis). Results for component e2 are similar and are not shown. Left-hand 
panels show individual simulated galaxies, right-hand panels show results 
binned in intervals of the measured ellipticity. Two magnitude ranges are 
shown, m > 22 (upper panels) and m ^ 22 (lower panels). The solid 
lines have a slope of unity, the dashed lines on the left-hand panels show 
the least-squares regression of input values on estimated values. The mean 
error on individual measured ellipticities is shown on the left-hand panels. 
Vertical error bars on the right-hand panels indicate the error in the mean 
input values in each interval of measured values. 

13 percent of galaxies in the STEP simulations. Some rejection of 
close pairs is required in most other methods of shape estimation 
also: in future, development of a fast multiple-component fitting al- 
gorithm might allow this constraint to be relaxed. We also test the 
fit returned by the fast fitting algorithm to determine whether the fit- 
ted centroid of a galaxy is within a reasonable range of the nominal 
position, given the prior on the galaxy position: this would iden- 
tify some of the cases of multiple galaxies. The criterion adopted is 
that the fitted galaxy position should lie within 3<r of the nominal 
position, where a 2 is the prior position variance, and this excludes 
9 percent of the initial simulated galaxy sample but no others that 
are excluded by the 'close pairs' criterion. 

Fig.[2]shows the posterior probability surfaces that result from 
fitting to two of the simulated galaxies. For completeness we also 
show the likelihood surfaces, which are broader and more biased 
away from the nominal value of ellipticity. 

Fig.[3] shows the results for each galaxy in the simulations 
(only the first component of ellipticity, ei, is shown, similar results 
are obtained for 62). At bright magnitudes there is good correspon- 
dence between input and measured ellipticity values. The slope ap- 
pears slightly steeper than unity, but with a value for the slope of 
1.04 ± 0.08 the departure from unity is not very significant. 

At fainter magnitudes, as the signal-to-noise decreases, an in- 
creasing fraction of galaxies with a given value of the Bayesian 
measure are drawn from a wider range of input ellipticities, as ex- 
pected from the earlier discussion. The slope of the relation be- 
tween input and measured values is again close to unity, with value 
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1.0 s per galaxy, implying that a survey of 10 7 galaxies could be 
analysed in a few months on a single standard desktop PC. The 
computation time scales inversely with the square of the sampling 
in ellipticity and increases approximately as m 2 log m for subim- 
ages of size m x m. 





5 FURTHER CONSIDERATIONS 
5.1 The ellipticity prior 

A number of studies have been made of the distribution 



Figure 4. The summed posterior probability distribution of measured ellip- 
ticity values ei (top) and e2 (bottom) as a function of apparent magnitude. 
The prior Vie) is also shown for comparison as a dashed line. The magni- 
tude ranges of the simulated galaxies are m $C 22 (left panels) and m > 22 
(right panels). 



0.93 ± 0.11. At all magnitudes the summed posterior probability 
distribution is a faithful reproduction of the distribution of the in- 
put prior distribution (Fig.Q as expected from section l272l There is 
also no detectable correlation between estimated values of ei and 
e2 in this simulated galaxy sample. 

We can also investigate the effect of measurement uncertainty 
in the prior positions of the galaxies. A random position uncertainty 
drawn from a normal distribution was introduced to each galaxy 
and the shapes remeasured. The prior assumed in the fitting was a 
normal distribution of rms 3 pixels throughout. In the results, no 
change was found in the slope of Fig.[3]for rms position uncertain- 
ties as large as 10 pixels. The scatter about the mean relation did 
not change for rms uncertainties less than 3 pixels and increased by 
4 percent for rms uncertainties as large as 10 pixels. This test indi- 
cates that the results are not sensitive to uncertainties in galaxy po- 
sition measurement. We would recommend that the position prior 
that is chosen should match the actual position uncertainty for the 
faintest galaxies that are reliably fitted. 

Tests of the algorithm, again using the full suite of STEP simu- 
lations specifically to measure the shear values recovered, are made 
in the companion paper (Kitching et al. in preparation). 



4.2 Speed 

The algorithm has been implemented in the C programming lan- 
guage □ for use on desktop computing systems, with discrete 
fast Fourier transform s being supplied by the FFTW library 
dFrigo & Stevensir2 005 ) B The computational speed per galaxy ob- 
viously depends on the computing system being used as well as 
on issues such as the extent to which the multiplex advantage of 
having many galaxies per PSF function can be exploited. In the 
simulations described above, using readily available 2 GHz desk- 
top PCs in 2007 and evaluating the likelihood surfaces on a grid 
of sampling interval 0.1 in e, we found computation times around 



2 The code lensfit is available on request from the authors: modification to 
the data input stages is likely to be required for any particular survey. 

3 http://www.fftw.org 



of galaxy ellipticities (e.g. lLam bas. Maddox & Lovedav 1992; 
Bramerd. Blandford & Smail 1 19961 : lEbbels. Kneib & Ellisl ll99S : 
Bernstein & Jarvisl 20021) . These studies find a wide variation in 
distribution of axial ratios, which appears strongly dependent on 
apparent magnitude, presumably largely as a result of the chang- 
ing mix of galaxies with brightness and redshift. The distribution 
of ellipticities at the faint magnitudes probed by ongoing and fu- 
ture weak lensing surveys is even less well-known, and the best 
estimate would come from the lensing data itself. For a sufficiently 
large survey the prior estimate could also be allowed to be a func- 
tion of galaxy brightness, redshift or colour, if that information 
were available. One way to estimate the ellipticity prior may be to 
adopt an iterative approach: evaluate the summed posterior proba- 
bility distribution starting from an initial guess of the prior distri- 
bution; then iteratively adjust the assumed prior until the summed 
posterior and prior distributions agree. We would expect this to be 
a stable iteration in the absence of sampling noise, because if the 
prior is initially assumed distributed to values that are smaller than 
are required to explain the data, the next iteration will adjust the 
prior to be distributed to large values, and vice versa. Such an ap- 
proach might however be unstable with small surveys where sam- 
pling noise might be important. 

In the case of lensing shear estimation, the ellipticity prior 
should also include the shear effect, and should not just be the in- 
trinsic pre-sheared distribution. As the shear varies on relatively 
small scales, and we are unlikely to have sufficient number of 
galaxies to measure accurately the ellipticity distribution in small 
regions, we suggest that correct generation of the prior should be to 
force the prior to be circularly symmetric, centred on (e) = 0, and 
to be obtained from the large numbers of galaxies that comprise the 
full survey. In this way 'false' shear variation arising from noise on 
the prior would be avoided, but the resulting shear values would be 
slightly biased to low values, in a magnitude-dependent way. This 
bias has already been discussed in section |2~4l and a method of cor- 
recting for the bias using the shear sensitivity has been described. 



5.2 Choice of model surface brightness profile 

A key advantage of the model-fitting approach over other methods 
is that a surface brightness profile may be chosen that accurately 
represents the actual profiles of galaxies. Two obvious choices of 
profile are exponential or de Vaucouleurs. In fact, it is notoriously 
difficult to choose between these profiles when fitting to faint galax- 
ies, so we do not expect the accuracy of the weak lensing mea- 
surement to depend strongly on which of these profiles is chosen. 
Some greater freedom in profile could be allowed by adding the 
Sersic index as a free parameter, allowing exponential and de Vau- 
couleurs m odels to be treated as special cases of this generalised 
profile (e.g. Dun lop et al.ll2003l) , however it is unlikely that the ad- 
dition of an extra parameter can be justified on evidence grounds. 
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A similar consideration is that galaxies generally are composed of 
bulge and disc components, which when viewed at an inclined an- 
gle may present differing ellipticities: accurate modelling of galax- 
ies requires both components to be fitted, but again, for generic 
shape measurement of faint galaxies where information at this level 
of sophistication is not present in the data, this seems unwarranted. 



5.3 Addition of multiple images or wavebands 

It may be that a weak lensing survey comprises multiple images of 
the same region of sky, but taken at different times and hence with 
differing PSFs, and possibly in different wavebands. The latter is 
likely if broad-band optical photometric redshifts are also being 
estimated from the same data that is being used for shape measure- 
ment. Clearly one would like to optimally estimate galaxy shape 
from the combination of all this data, but it would not be optimal 
simply to co-add the data, because of the differing PSFs in each im- 
age. The model-fitting algorithm described above allows a natural 
way to optimally include all the data, since all we need do is add 
the likelihoods for the models fitted to each galaxy. In doing this, 
we should take care that the nominal galaxy positions are the same 
in each image, so the optimal way to proceed would be to co-add 
images for the purpose of detecting and measuring nominal galaxy 
positions only, and then fitting each individual image with models 
convolved with the appropriate PSF and adding the resulting likeli- 
hoods. Images with a mixture of seeing qualities are thus optimally 
combined for the shape measurement. 



5.4 Weak lensing from radio interferometer data 

It is clear that in large future optical surveys systematic uncertain- 
ties in PSF correction will be a dominating concern, indeed this is 
a significant factor in the case for space-based weak lensing mis- 
sions. Ground-based optical PSFs vary temporally and very often 
on spatial scales comparable to those on which the cosmic shear 
signal is detectable. Even HST lensing s t udies suffer significantl y 
from PSF variation ^Rhodes et afll2007t ISchrabback et alj|2007t) . 
In principle radio interferometers have precisely known PSFs, be- 
ing determined by the antenna positions (note that full 3D knowl- 
edge of antenna positions is required, to allow for curvature of the 
Earth and natural height variations). The PSF varies with hour an- 
gle and declination, but in a completely deterministic way. Other 
effects such as bandwidth and sampling-time smearing can also be 
precisel y computed and inc orporated into the shape measurement 
process dChang et alj|2004h . Because interferometer measurement 
are made in the Fourier domain, and because the noise also orig- 
inates in that domain (being associated with individual antennas) 
it makes sense to measure galaxy shapes in Fourier space (in the 
image plane the noise is correlated between pixels, eff e ctivel y be- 
ing also convolved with the PSF). IChang & Refregier (2002) and 
Chang et al.l d2004[> have already shown how a shapelets dRefregier 
2003URefregier & BacorJ [2003) based approach can be extended 
to the Fourier domain. The Bayesian algorithm presented in this 
paper already operates in the Fourier domain, so it should be eas- 
ily adapted for radio interferometer data, which will be particularly 
relevant for future deep radio surveys such as those proposed for 
the Square Kilometer Array. 



6 CONCLUSIONS 

We have argued that a model-fitting approach to galaxy shape mea- 
surement should provide an optimum approach to shape measure- 
ment for large weak-lensing surveys, with the advantages that the 
signal-to-noise of the shape measurement should be optimised and 
random measurement errors can be estimated. We have further ar- 
gued that a Bayesian estimation process allows unbiased shape es- 
timation to be made, although even in a realistic implementation of 
a Bayesian method there is a bias in recovered shear values intro- 
duced by the presence of the prior probability distribution. This bias 
may be calculated from the measured likelihood surfaces, however, 
and in this paper we have spent some time discussing the calcula- 
tion of the 'shear sensitivity'. Overall this approach to shape mea- 
surement should provide a framework for shear measurement that 
does not need external calibration by comparison with simulations. 

A traditional disadvantage of model-fitting is that it may be 
computationally time-consuming, and in this paper we present a 
fast algorithm for measuring the shapes of individual galaxies. 
The algorithm makes use of analytic marginalisation over surface 
brightness amplitude, and by working in Fourier space enables 
rapid marginalisation over galaxy position. The algorithm has been 
tested and has an adequate speed on current generations of comput- 
ers for use with large ongoing and planned weak-lensing surveys. 
Close pairs of galaxies are not treated by the algorithm, but pro- 
vided such close pairs can be identified in the data a separate fitting 
process may be applied to those. 

The Bayesian method and fast fitting algorithm have been 
tested o n simulated galaxies created for the Shear T Esting Program 
(STEP: lHevmans etal. 20061: iMassev et al.| [2007b) and promising 
results on the measurement of individual galaxy ellipticities have 
been obtained. A companion paper (Kitching et al. in preparation) 
will test the measurement of weak lensing shear in the STEP simu- 
lations. 
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